Application for chapter 4 of the thesis
Covariates impacts in SAR models for CoDa
C. Thomas-Agnan, T. Laurent and A. Ruiz-Gazen
Last update: 2025-09-22
This document provides the R code used to reproduce the results presented in Chapter 4 of Thibault Laurent’s thesis. The methodological framework is based on the original article, but the empirical application has been adapted here to address the migration–climate context developed in the dissertation. The original codes corresponding to the published paper are available at: https://github.com/tibo31/spatial_coda_elasticities
Packages needed:
library(cartography) # mapping
library(compositions) # compositions data
library(gridExtra) # several frame in ggplot
library(Matrix) # sparse matrix
library(scatterpie) # plot several pie
library(sp) # spatial package sp
library(spdep) # spatial econometrics
library(sf) # spatial package sf
library(tidyverse) # tidyverse universe Import data
Flow data
We import first the flow data estimated in Chapter 1, using the optimal transport approach, which serve as the dependent variable in the spatial compositional regression models of Chapter 4.
load(file = "data/our_estimates_2024.RData")
temp <- our_estimates_5y %>%
filter(type == "tot",
year == "2020-2024") %>%
select(origin, dest, hat_transport_open_outwards, hat_transport_open_returns, hat_transport_open_transits)
outflows <- temp %>%
group_by(origin) %>%
summarize(out_outwards = sum(hat_transport_open_outwards),
out_returns = sum(hat_transport_open_returns),
out_transit = sum(hat_transport_open_transits)) %>%
mutate(out_total = out_outwards + out_returns + out_transit) %>%
rename(iso3 = origin)
inflows <- temp %>%
group_by(dest) %>%
summarize(in_outwards = sum(hat_transport_open_outwards),
in_returns = sum(hat_transport_open_returns),
in_transit = sum(hat_transport_open_transits)) %>%
mutate(in_total = in_outwards + in_returns + in_transit) %>%
rename(iso3 = dest)Spatial contours
We import the spatial contours of countries, which will be used to match the flow data with geographical units.
my_proj <- '+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
source(file = "data/Load_Geospatial_Data.R")
sf_use_s2(F)## Spherical geometry (s2) switched offworld_sf <- GeoDATA %>%
rename(region = SIREGION, iso3 = ISO3)
countries_regions <- world_sf %>%
rmapshaper::ms_simplify(, keep = 0.05, keep_shapes = T)
# polygons with same country than UN data basis
world <- st_transform(countries_regions, my_proj)
world_union <- st_union(world)
# contours, longitude and latitude of the map
p1 <- rbind(c(-180, -90),
c(180, -90),
cbind(180, seq(-90, 90, by = 1)),
c(180, 90),
c(-180, 90),
cbind(-180, seq(90, -90, by = -1)),
c(-180, -90))
border_sf <- st_sfc(st_polygon(list(p1)))
long_sf <- st_sfc(st_linestring(cbind(-150, seq(-90, 90, by = 1))),
st_linestring(cbind(-120, seq(-90, 90, by = 1))),
st_linestring(cbind(-90, seq(-90, 90, by = 1))),
st_linestring(cbind(-60, seq(-90, 90, by = 1))),
st_linestring(cbind(-30, seq(-90, 90, by = 1))),
st_linestring(cbind(0, seq(-90, 90, by = 1))),
st_linestring(cbind(150, seq(-90, 90, by = 1))),
st_linestring(cbind(120, seq(-90, 90, by = 1))),
st_linestring(cbind(90, seq(-90, 90, by = 1))),
st_linestring(cbind(60, seq(-90, 90, by = 1))),
st_linestring(cbind(30, seq(-90, 90, by = 1)))
)
lat_sf <- st_sfc(st_linestring(cbind(seq(-180, 180, by = 1), -60)),
st_linestring(cbind(seq(-180, 180, by = 1), -30)),
st_linestring(cbind(seq(-180, 180, by = 1), -0)),
st_linestring(cbind(seq(-180, 180, by = 1), 30)),
st_linestring(cbind(seq(-180, 180, by = 1), 60))
)
# Coordinate Reference System
st_crs(border_sf) <- 4326
st_crs(long_sf) <- 4326
st_crs(lat_sf) <- 4326
border_sf <- st_transform(border_sf, my_proj)
long_sf <- st_transform(long_sf, my_proj)
lat_sf <- st_transform(lat_sf, my_proj)
sea <- st_difference(border_sf, world_union)plot(sea, col = "lightblue", border = rgb(0.2, 0.2, 0.2))
plot(long_sf, add = T, col = "lightgrey")
plot(lat_sf, add = T, col = "lightgrey")
We create the object that will gather all variables,
combining the covariates imported previously with the new stock-based
indicators.
world_iso_3 <- merge(world, outflows, by = "iso3")
world_iso_3 <- merge(world_iso_3, inflows, by = "iso3")Explanatory variables
We import the same set of covariates as in Chapter 3,
so that the analysis remains comparable across chapters.
load("data/covariates.RData")
world_iso_3 <- merge(world_iso_3, my_X[, -c(2:37)], by = "iso3")
world_iso_3$class_conflict <- factor(world_iso_3$class_conflict,
levels = c("low", "medium", "high", "very high"))We add a new variable based on migration stocks,
namely the percentage of migrants residing in the country
and the percentage of nationals living abroad.
load("data/stock_tot_2020.RData")
stock_info <- data.frame(
iso3 = colnames(s1),
percent_foreign = 1 - diag(s1_open) / colSums(s1_open),
percent_outside = 1 - diag(s1_open) / rowSums(s1_open)
)
world_iso_3 <- merge(world_iso_3, stock_info, by = "iso3")Exploratory analysis
Dependent variable
Tricolore map
We begin by exploring the dependent variable, which corresponds to the compositional structure of migration flows (outward, return, transit). Here, we reproduce Figure 4.6 of the article: a ternary diagram with colors and mapping, which provides a first visualization of the distribution of flows across the three components.
library(compositions)
A <- c(0, 0)
B <- c(1, 0)
C <- c(0.5, sqrt(3) / 2)
# identify countries influenced by one composition
# out
ind_belize <- 25
ind_saudi <- 177
ind_viet <- which(world_iso_3$SOVRN == "Viet Nam")
temp_out <- world_iso_3[c(ind_belize, ind_saudi, ind_viet), ]
ind_chad <- which(world_iso_3$SOVRN == "Chad")
ind_cuba <- which(world_iso_3$SOVRN == "Cuba")
ind_lebanon <- which(world_iso_3$SOVRN == "Lebanon")
temp_in <- world_iso_3[c(ind_lebanon, ind_cuba, ind_chad), ]
temp_compo <- st_drop_geometry(world_iso_3)
# merge spatial data
for(s in c("out", "in")) {
s_3 <- temp_compo[, paste0(s, "_", c("outwards", "returns", "transit"))] /
temp_compo[, paste0(s, "_total")]
s_3[s_3 == 0] <- 0.01
acomp_s3 <- compositions::acomp(s_3)
names(acomp_s3) <- c("outwards", "returns", "transit")
n <- nrow(s_3)
s_3_mean <- acomp_s3[1, ] / n
for(k in 2:nrow(acomp_s3))
s_3_mean <- s_3_mean + acomp_s3[k, ] / n
s_3_centered <- as(acomp_s3, "matrix")
Y_simplex_x <- s_3_centered[, 1] * A[1] + s_3_centered[, 2] * B[1] +
s_3_centered[, 3] * C[1]
Y_simplex_y <- s_3_centered[, 1] * A[2] + s_3_centered[, 2] * B[2] +
s_3_centered[, 3] * C[2]
#pdf(paste0("figures/coda_colors_", s , ".pdf"),
# width = 13 * 0.8, height = 4.35 * 0.8)
par(oma = c(0, 0, 0, 0), mfrow = c(1, 2))
par(mar = c(0, 0, 0.3, 0))
pal_col_ternary <- tricolore::Tricolore(data.frame(
outwards = c(0.7, 0.6, 0.4, 0.25, 0.05, 0.5, 0.25, 0.05, 0.25),
returns = c(0.25, 0.25, 0.25, 0.25, 0.25, 0.45, 0.45, 0.45, 0.7),
transit = c(0.05, 0.15, 0.35, 0.5, 0.7, 0.05, 0.2, 0.5, 0.05)),
'outwards', 'returns', 'transit',
hue = 0.2, breaks = 3)$rgb
list_poly <- list(
list(x = c(0, 1/3, 1/6, 0),
y = c(0, 0, C[2] / 3, 0)),
list(x = c(1/6, 1/3, 1/2, 1/6),
y = c(C[2]/3, 0, C[2]/3, C[2]/3)),
list(x = c(1/3, 2/3, 1/2, 1/3),
y = c(0, 0, C[2]/3, 0)),
list(x = c(1/2, 2/3, 5/6, 1/2),
y = c(C[2]/3, 0, C[2]/3, C[2]/3)),
list(x = c(2/3, 1, 5/6, 2/3),
y = c(0, 0, C[2]/3, 0)),
list(x = c(1/6, 1/2, 1/3, 1/6),
y = c(C[2]/3, C[2]/3, 2*C[2]/3, C[2]/3)),
list(x = c(1/3, 1/2, 2/3, 1/3),
y = c(2*C[2]/3, C[2]/3, 2*C[2]/3, 2*C[2]/3)),
list(x = c(1/2, 5/6, 2/3, 1/2),
y = c(C[2]/3, C[2]/3, 2 * C[2]/3, C[2]/3)),
list(x = c(1/3, 2/3, 1/2, 1/3),
y = c(2 * C[2]/3, 2 * C[2]/3, C[2], 2 * C[2]/3))
)
plot(rbind(A, B, C), xaxt = "n", yaxt = "n", frame = F,
type = "n", xlab = "", ylab = "", asp = 1, main = "",
ylim = c(-0.12, sqrt(3)/2), cex.main = 2)
for(l in 1:9) {
polygon(list_poly[[l]]$x, list_poly[[l]]$y, col = pal_col_ternary[l],
border = NULL, lty = 1, lwd = 0.1)
}
points(Y_simplex_x, Y_simplex_y, col = "black",
cex = 0.4, pch = 16)
if(s == "out") {
points(Y_simplex_x[c(ind_belize, ind_saudi, ind_viet)],
Y_simplex_y[c(ind_belize, ind_saudi, ind_viet)],
col = "red",
cex = 0.6, pch = 16)
} else {
points(Y_simplex_x[c(ind_lebanon, ind_cuba, ind_chad)],
Y_simplex_y[c(ind_lebanon, ind_cuba, ind_chad)],
col = "red",
cex = 0.6, pch = 16)
}
text(c(0.04, 0.96, 1/2 - 0.1), c(0-0.12, -0.12, sqrt(3)/2 - 0.041),
c("Outward", "Return", "Transit"), pos = 3,
cex = 1, col = pal_col_ternary[c(1, 5, 9)])
lines(c(0, 1), c(0, 0), col = "black")
lines(c(0, 0.5), c(0, sqrt(3)/2), col = "black")
lines(c(1, 0.5), c(0, sqrt(3)/2), col = "black")
mtext(paste0("\t \t \t Compositions of ", s, "flows"),
side = 2, line = -2, cex = 1.1)
# XR lines
my_labels <- c("1/3", "2/3")
for (k in 1:2) {
text(1-k/6, k/3 * C[2], my_labels[k], pos = 4, cex = 0.7,
col = pal_col_ternary[c(9)])
}
# Left lines
for (k in 1:2) {
text(k/3, 0, my_labels[k], pos = 1, cex = 0.7, srt = 250,
col = pal_col_ternary[c(5)])
}
# Right lines
for (k in 1:2) {
text(k * 1 / 6-0.01, k/3 * C[2], rev(my_labels)[k], pos = 3, cex = 0.7, srt = 135,
col = pal_col_ternary[c(1)])
}
new_col <- tricolore::Tricolore(as(acomp_s3, "data.frame"),
'outwards', 'transit','returns',
hue = 0.2, breaks = 3)$rgb
#####################################.
par(mar = c(0, 0, 0, 0))
plot(sea, col = "lightblue", border = rgb(0.2, 0.2, 0.2))
plot(long_sf, add = T, col = "lightgrey")
plot(lat_sf, add = T, col = "lightgrey")
plot(st_geometry(world_iso_3), border = new_col,
lwd = 0.2,
col = new_col, add = T)
if(s == "out") {
plot(st_geometry(temp_out),
border = "red",
col = pal_col_ternary[c(9, 5, 1)],
lwd = 0.5, add = T)
text(st_coordinates(st_centroid(temp_out))[, 1],
st_coordinates(st_centroid(temp_out))[, 2],
temp_out$SOVRN, pos = 4, cex = 0.5)
} else {
plot(st_geometry(temp_in),
border = "red",
col = pal_col_ternary[c(9, 5, 1)],
lwd = 0.5, add = T)
text(st_coordinates(st_centroid(temp_in))[, 1],
st_coordinates(st_centroid(temp_in))[, 2],
temp_in$SOVRN, pos = 4, cex = 0.5)
}
#dev.off()
}## Registered S3 methods overwritten by 'ggtern':
## method from
## grid.draw.ggplot ggplot2
## plot.ggplot ggplot2
## print.ggplot ggplot2## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
## Warning: st_centroid assumes attributes are constant over geometries
We now examine the countries with the highest shares of outward migration in their compositional structure. This allows us to identify which origins contribute most strongly to departure flows.
temp_compo %>%
mutate(out_outwards = out_outwards / out_total,
out_returns = out_returns / out_total,
out_transit = out_transit / out_total,
in_outwards = in_outwards / in_total,
in_returns = in_returns / in_total,
in_transit = in_transit / in_total) %>%
arrange(-out_outwards) %>%
head(30)Barplot
As an alternative to the ternary diagram, we represent the
compositions using barplots
(see Figure D.2 in the Appendix). This view highlights more clearly the
relative weight of outward, return, and transit flows for each
country.
# inflows
temp_compo %>%
mutate(in_outward = in_outwards / in_total,
in_return = in_returns / in_total,
in_transit = in_transit / in_total) %>%
mutate(iso3 = fct_reorder(iso3, in_outward)) %>%
select(iso3, CONTINENT, in_outward, in_return, in_transit) %>%
pivot_longer(
cols = c("in_outward", "in_return", "in_transit"),
names_to = "type",
values_to = "share") %>%
mutate(type = substr(type, 4, nchar(type))) %>%
ggplot() +
geom_col(aes(x = iso3, fill = type, y = share)) +
theme(axis.text.x = element_text(angle = 90)) +
xlab("Country") + ylab("Shares of migrants in inflows") +
facet_grid(. ~ CONTINENT, scales = "free_x", space = "free", ) +
scale_fill_grey()
ggsave("figures/barplot_in.pdf", width = 25, height = 5)
# outflows
temp_compo %>%
mutate(out_outward = out_outwards / out_total,
out_return = out_returns / out_total,
out_transit = out_transit / out_total) %>%
mutate(iso3 = fct_reorder(iso3, out_outward)) %>%
select(iso3, CONTINENT, out_outward, out_return, out_transit) %>%
pivot_longer(
cols = c("out_outward", "out_return", "out_transit"),
names_to = "type",
values_to = "share") %>%
mutate(type = substr(type, 4, nchar(type))) %>%
ggplot() +
geom_col(aes(x = iso3, fill = type, y = share)) +
theme(axis.text.x = element_text(angle = 90)) +
xlab("Country") + ylab("Shares of migrants in outflows") +
facet_grid(. ~ CONTINENT, scales = "free_x", space = "free", ) +
scale_fill_grey()
ggsave("figures/barplot_out.pdf", width = 25, height = 5)Spatial weight matrix
We represent the spatial weight matrix used in the analysis, which captures the neighborhood structure between countries.
world_iso_3_agg <- world_iso_3
world_iso_3_agg$iso3[world_iso_3_agg$iso3 == "ESH"] <- "MAR"
world_iso_3_agg <- aggregate(world_iso_3_agg[, "iso3"],
by = list(iso3 = world_iso_3_agg$iso3),
FUN = length)
# fontiere Maroc / Sahara
frontiers <- st_intersection(world_iso_3[world_iso_3$iso3 == "MAR", ],
world_iso_3[world_iso_3$iso3 == "ESH", ])my_centroid <- st_as_sf(data.frame(
long = world_iso_3$Centroid_of_LARGEST_POLYGON_X,
lat = world_iso_3$Centroid_of_LARGEST_POLYGON_Y,
iso3 = world_iso_3$iso3),
coords = c("long", "lat"),
crs = 4326)
# based on contiguity + 3 neighbours
nb_out <- union.nb(poly2nb(world_iso_3),
knn2nb(knearneigh(my_centroid, k = 3)))
OW <- nb2mat(nb_out)
crs_robin <- st_crs(world_iso_3)
neighbours_sf <- st_as_sf(nb2lines(nb_out,
coords = st_coordinates(st_centroid(world_iso_3))))
st_crs(neighbours_sf) <- crs_robin
long_sf_180 <- st_sfc(st_linestring(cbind(-180, seq(-90, 90, by = 1))),
st_linestring(cbind(180, seq(-90, 90, by = 1))),
crs = 4326
)
long_sf_180 <- st_transform(long_sf_180, crs_robin)
ind_big <- which(as.numeric(st_length(neighbours_sf)) > 10^7)
#pdf("figures/spatial_weight.pdf", width = 10, height =5)
par(oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0))
plot(st_geometry(world_iso_3_agg), border = "lightgrey",
ylim = c(-7*10^6, 7*10^6))
plot(st_geometry(frontiers), add = T, col = "lightgrey",
lty = 3)
plot(st_geometry(sea), col = "lightblue", add = T)
plot(st_geometry(neighbours_sf[-ind_big, ]), col = "darkred",
add =T, lwd = 0.8)
for(l in 1:length(ind_big)) {
# point 1
x_pt_1 <- neighbours_sf[ind_big[l], ]$geometry[[1]][1, 1]
y_pt_1 <- neighbours_sf[ind_big[l], ]$geometry[[1]][1, 2]
pt_1 <- st_sfc(st_point(c(x_pt_1, y_pt_1)), crs = crs_robin)
# point 2
x_pt_2 <- neighbours_sf[ind_big[l], ]$geometry[[1]][2, 1]
y_pt_2 <- neighbours_sf[ind_big[l], ]$geometry[[1]][2, 2]
pt_2 <- st_sfc(st_point(c(x_pt_2, y_pt_2)), crs = crs_robin)
seg1 <- st_nearest_points(pt_1, long_sf_180[which.min(as.numeric(st_distance(pt_1, long_sf_180))), ])
seg2 <- st_nearest_points(pt_2, long_sf_180[which.min(as.numeric(st_distance(pt_2, long_sf_180))), ])
plot(st_geometry(seg1), add = T, col = "darkred")
plot(st_geometry(seg2), add = T, col = "darkred")
}
Scatter plots
We plot the scatter plots showing how the components of the dependent variable (outward, return, transit) relate to each of the explanatory covariates. This provides a first visual assessment of possible associations or patterns.
nom_vars <- c("hwf", "dry", "Lpop", "HDI", "politicalstability",
"percent_foreign", "percent_outside")
names_title <- c(
"Heat wave frequency",
"Consecutive dry days",
"Log population",
"Human Development Index",
"Political stability",
"Percent foreign residents",
"Percent natives living abroad"
)
# Quantitative
for(l in 1:length(nom_vars)) {
k <- nom_vars[l]
#pdf(paste0("figures/explo_", k, ".pdf"), width = 8, height = 5)
par(mfrow = c(2, 3), oma = c(0, 0, 0.2, 0), mar = c(3.3, 4, 1, 0.5), las = 1,
mgp = c(2.2, 1, 0))
for(s in c("out", "in")) {
plot(world_iso_3[[k]], world_iso_3[[paste0(s, "_outwards")]] /
world_iso_3[[paste0(s, "_total")]],
xlab = names_title[l], ylab = paste0("Share of outwards among ", s, "flows"),
pch = 16, col = pal_col_ternary[c(1)], ylim = c(0, 1))
abline(lm(world_iso_3[[paste0(s, "_outwards")]] / world_iso_3[[paste0(s, "_total")]] ~
world_iso_3[[k]]), col = "red")
if(s == "out")
title("Outward")
plot(world_iso_3[[k]], world_iso_3[[paste0(s, "_returns")]] / world_iso_3[[paste0(s, "_total")]],
xlab = names_title[l], ylab = paste0("Shares of returns among ", s, "flows"),
pch = 16, col = pal_col_ternary[c(5)], ylim = c(0, 1))
abline(lm(world_iso_3[[paste0(s, "_returns")]] / world_iso_3[[paste0(s, "_total")]] ~
world_iso_3[[k]]), col = "red")
if(s == "out")
title("Return")
plot(world_iso_3[[k]], world_iso_3[[paste0(s, "_transit")]] / world_iso_3[[paste0(s, "_total")]],
xlab = names_title[l], ylab = paste0("Shares of transits among ", s, "flows"),
pch = 16, col = pal_col_ternary[c(9)], ylim = c(0, 1))
abline(lm(world_iso_3[[paste0(s, "_transit")]] / world_iso_3[[paste0(s, "_total")]] ~
world_iso_3[[k]]), col = "red")
if(s == "out")
title("Transit")
}
#dev.off()
}
noms_out <- c("out_outwards", "out_returns", "out_transit",
"dry", "hwf", "politicalstability",
"HDI", "Lpop", "percent_foreign", "percent_outside")
noms_in <- c("in_outwards", "in_returns", "in_transit", "dry", "hwf",
"politicalstability",
"HDI", "Lpop", "percent_foreign", "percent_outside")
cor(st_drop_geometry(world_iso_3)[, noms_out])## out_outwards out_returns out_transit dry hwf
## out_outwards 1.00000000 0.18692550 0.20018428 0.06190276 -0.08711607
## out_returns 0.18692550 1.00000000 0.98968595 0.04286184 -0.14226295
## out_transit 0.20018428 0.98968595 1.00000000 0.05464424 -0.15733734
## dry 0.06190276 0.04286184 0.05464424 1.00000000 -0.01581172
## hwf -0.08711607 -0.14226295 -0.15733734 -0.01581172 1.00000000
## politicalstability -0.42790637 -0.06353307 -0.06218877 -0.40332414 0.20058276
## HDI -0.11686190 0.24332799 0.27041782 -0.24371845 0.04388493
## Lpop 0.42592193 0.36632032 0.39667132 0.22424889 -0.35274096
## percent_foreign -0.19960143 0.06749276 0.08026922 0.08442611 0.31015859
## percent_outside -0.05671371 -0.16117649 -0.17552744 -0.14068704 0.17387081
## politicalstability HDI Lpop percent_foreign
## out_outwards -0.42790637 -0.11686190 0.4259219 -0.19960143
## out_returns -0.06353307 0.24332799 0.3663203 0.06749276
## out_transit -0.06218877 0.27041782 0.3966713 0.08026922
## dry -0.40332414 -0.24371845 0.2242489 0.08442611
## hwf 0.20058276 0.04388493 -0.3527410 0.31015859
## politicalstability 1.00000000 0.63885426 -0.6278686 0.45692080
## HDI 0.63885426 1.00000000 -0.2209145 0.50374807
## Lpop -0.62786858 -0.22091449 1.0000000 -0.47011242
## percent_foreign 0.45692080 0.50374807 -0.4701124 1.00000000
## percent_outside 0.17367858 0.08819054 -0.4643476 0.13696442
## percent_outside
## out_outwards -0.05671371
## out_returns -0.16117649
## out_transit -0.17552744
## dry -0.14068704
## hwf 0.17387081
## politicalstability 0.17367858
## HDI 0.08819054
## Lpop -0.46434764
## percent_foreign 0.13696442
## percent_outside 1.00000000cor(st_drop_geometry(world_iso_3)[, noms_in])## in_outwards in_returns in_transit dry hwf
## in_outwards 1.00000000 0.28443860 0.95587598 0.04040555 -0.15096361
## in_returns 0.28443860 1.00000000 0.28248827 0.09984444 -0.16981585
## in_transit 0.95587598 0.28248827 1.00000000 0.04848808 -0.17003457
## dry 0.04040555 0.09984444 0.04848808 1.00000000 -0.01581172
## hwf -0.15096361 -0.16981585 -0.17003457 -0.01581172 1.00000000
## politicalstability -0.04983782 -0.47588659 -0.03212315 -0.40332414 0.20058276
## HDI 0.23941252 -0.09816766 0.29046146 -0.24371845 0.04388493
## Lpop 0.36630221 0.56441731 0.38854119 0.22424889 -0.35274096
## percent_foreign 0.04597390 -0.27253176 0.08499821 0.08442611 0.31015859
## percent_outside -0.16416912 -0.04441739 -0.17488521 -0.14068704 0.17387081
## politicalstability HDI Lpop percent_foreign
## in_outwards -0.04983782 0.23941252 0.3663022 0.04597390
## in_returns -0.47588659 -0.09816766 0.5644173 -0.27253176
## in_transit -0.03212315 0.29046146 0.3885412 0.08499821
## dry -0.40332414 -0.24371845 0.2242489 0.08442611
## hwf 0.20058276 0.04388493 -0.3527410 0.31015859
## politicalstability 1.00000000 0.63885426 -0.6278686 0.45692080
## HDI 0.63885426 1.00000000 -0.2209145 0.50374807
## Lpop -0.62786858 -0.22091449 1.0000000 -0.47011242
## percent_foreign 0.45692080 0.50374807 -0.4701124 1.00000000
## percent_outside 0.17367858 0.08819054 -0.4643476 0.13696442
## percent_outside
## in_outwards -0.16416912
## in_returns -0.04441739
## in_transit -0.17488521
## dry -0.14068704
## hwf 0.17387081
## politicalstability 0.17367858
## HDI 0.08819054
## Lpop -0.46434764
## percent_foreign 0.13696442
## percent_outside 1.00000000We use boxplots to represent how the compositions of the dependent variable (outward, return, transit) vary across categories of the qualitative covariates.
#pdf("figures/boxplot.pdf", width = 10, height = 6)
my_conflict <- world_iso_3$class_conflict
par(mfrow = c(2, 3), oma = c(0, 0, 0, 0), mar = c(3, 4, 1, 1), las = 1)
for(s in c("out", "in")) {
# qualitative
boxplot(world_iso_3[[paste0(s, "_outwards")]] / world_iso_3[[paste0(s, "_total")]] ~
my_conflict, xlab = "conflict", ylab = paste0("Shares of ", s, "flows"),
main = "Outward")
boxplot(world_iso_3[[paste0(s, "_returns")]] / world_iso_3[[paste0(s, "_total")]] ~
my_conflict, xlab = "conflict", ylab = paste0("Shares of ", s, "flows"),
main = "Return")
boxplot(world_iso_3[[paste0(s, "_transit")]] / world_iso_3[[paste0(s, "_total")]] ~
my_conflict, xlab = "conflict", ylab = paste0("Shares of ", s, "flows"),
main = "Transit")
}
# dev.off()Modelling
OLS
We first run a simple linear regression on each component of the
ILR-transformed compositions.
This step serves as a baseline model, allowing us to assess the main
effects of the covariates
before introducing spatial dependence structures.
ilr_base <- ilrBase(D = 3)
v1 <- ilr_base[, 1]
v2 <- ilr_base[, 2]
ilr_base[1, 1] <- v2[3]
ilr_base[2:3, 1] <- v2[1:2]
ilr_base[1, 2] <- 0
ilr_base[2:3, 2] <- v1[1:2]
var_controls <- c("Lpop", "HDI", "politicalstability", "class_conflict",
"percent_foreign", "percent_outside")
vars_included <- c("hwf", "dry")
var_climate <- paste0(c(""), rep(vars_included, each = 1))
temp_lm <- st_drop_geometry(world_iso_3)
my_lm <- vector("list", 2)
names(my_lm) <- c("out", "in")
my_lm[["out"]] <- my_lm[["in"]] <- vector("list", 2)
names(my_lm[["out"]]) <- names(my_lm[["in"]]) <- c("ilr_1", "ilr_2")
for(type in c("out", "in")) {
temp_ilr <- ilr(temp_lm[, paste0(type, "_", c("outwards", "returns", "transit"))],
V = ilr_base)
temp_lm[[paste0(type, "_ilr_1")]] <- temp_ilr[, 1]
temp_lm[[paste0(type, "_ilr_2")]] <- temp_ilr[, 2]
form_str <- paste(paste0(type, "_ilr_1 ~"), paste(c(var_climate, var_controls),
collapse = " + "))
form_str_2 <- paste(paste0(type, "_ilr_2 ~"), paste(c(var_climate, var_controls),
collapse = " + "))
my_lm[[type]][["ilr_1"]] <- lm(as.formula(form_str), data = temp_lm)
my_lm[[type]][["ilr_2"]] <- lm(as.formula(form_str_2), data = temp_lm)
}Computing the predictions
We now compute the model predictions and visualize them. This step reproduces Figures D.10 and D.11, which compare observed and predicted values for the different components of the migration flow compositions.
new_data <- lapply(temp_lm[, c(var_climate, var_controls)],
function(x) ifelse(is.numeric(x), median(x), names(which.max(table(x))))) |>
as.data.frame()
for(type in c("out", "in")) {
#pdf(paste0("figures/predicted_OLS_", type, "flows.pdf"), width = 10, height = 5)
par(mfrow = c(2, 4), oma = c(0, 0, 0, 0), mar = c(3.5, 4, 0.5, 0.5), las = 1,
mgp = c(2.2, 1, 0))
for(s in 1:length(nom_vars)){
k <- nom_vars[s]
my_seq <- seq(min(temp_lm[[k]]), max(temp_lm[[k]]), length.out = 100)
res_pred <- matrix(0, 100, 2)
for(l in 1:100) {
new_data[[k]] <- my_seq[l]
predict_lm_1 <- predict(my_lm[[type]][["ilr_1"]], newdata = new_data)
predict_lm_2 <- predict(my_lm[[type]][["ilr_2"]], newdata = new_data)
res_pred[l, ] <- c(predict_lm_1, predict_lm_2)
}
res_pred_y <- ilrInv(res_pred, V = ilr_base)
pal_col <- RColorBrewer::brewer.pal(3, "Set2")
plot(my_seq, res_pred_y[, 1], type = "l", ylim = c(0, 1),
col = pal_col[1], ylab = paste0("Predicted shares in ", type, "flows"),
xlab = names_title[s])
lines(my_seq, res_pred_y[, 2], col = pal_col[2])
lines(my_seq, res_pred_y[, 3], col = pal_col[3])
if(s == 1)
legend("topleft", legend = c("Outward", "Return", "Transit"), lty = 1,
col = pal_col)
}
#dev.off()
}
Compute the semi-elasticities
We compute the semi-elasticities associated with the covariates, which quantify the relative change in each component of the migration composition resulting from a marginal variation of the explanatory variables.
res_impacts <- data.frame(
var = character(0),
model = character(0),
iso3 = character(0),
se1 = numeric(0),
se2 = numeric(0),
se3 = numeric(0))
add_unit <- c(hwf = 1, dry = 1, HDI = 0.00593,
politicalstability = 0.0443764,
percent_foreign = 0.01, percent_outside = 0.01)
for(type in c("in", "out")) {
for(k in nom_vars) {
# add unit
my_X_plus_only_ind_cnty <- st_drop_geometry(world_iso_3)
if(k != "Lpop") {
my_X_plus_only_ind_cnty[, k] <- my_X_plus_only_ind_cnty[, k] + add_unit[k]
} else {
# elasticities for Lpop
my_X_plus_only_ind_cnty[, "Lpop"] <- my_X_plus_only_ind_cnty[, "Lpop"] * 1.01
}
predict_lm_1 <- predict(my_lm[[type]][["ilr_1"]],
newdata = st_drop_geometry(world_iso_3))
predict_lm_2 <- predict(my_lm[[type]][["ilr_2"]],
newdata = st_drop_geometry(world_iso_3))
Y_pred <- ilrInv(cbind(predict_lm_1, predict_lm_2), V = ilr_base)
hat_y_ilr1 <- predict(my_lm[[type]][["ilr_1"]],
newdata = my_X_plus_only_ind_cnty)
hat_y_ilr2 <- predict(my_lm[[type]][["ilr_2"]],
newdata = my_X_plus_only_ind_cnty)
Y_plus <- ilrInv(cbind(hat_y_ilr1, hat_y_ilr2), V = ilr_base)
temp_diff <- as(Y_plus, "matrix") - as(Y_pred, "matrix")
names(temp_diff) <- paste0("se", 1:3)
res_impacts <- rbind(res_impacts,
data.frame(var = k,
model = type,
iso3 = world_iso_3$iso3,
se1 = temp_diff[, 1],
se2 = temp_diff[, 2],
se3 = temp_diff[, 3]))
}
}
res_impacts %>%
filter(model == "out",
iso3 %in% c("BLZ", "SAU", "VNM"))## var model iso3 se1 se2 se3
## 257 hwf out BLZ 0.0012852699 -0.0008866460 -3.986239e-04
## 1777 hwf out SAU 0.0003862733 -0.0004093371 2.306374e-05
## 2237 hwf out VNM 0.0010388637 -0.0007321584 -3.067052e-04
## 258 dry out BLZ -0.0005883713 0.0004048131 1.835582e-04
## 1778 dry out SAU -0.0001763973 0.0001856410 -9.243682e-06
## 2238 dry out VNM -0.0004773938 0.0003361330 1.412608e-04
## 259 Lpop out BLZ 0.0043812315 -0.0018080906 -2.573141e-03
## 1779 Lpop out SAU 0.0017168551 0.0001402927 -1.857148e-03
## 2239 Lpop out VNM 0.0048778472 -0.0029319296 -1.945918e-03
## 2510 HDI out BLZ -0.0048426608 0.0033620830 1.480578e-03
## 17710 HDI out SAU -0.0014456105 0.0015573830 -1.117725e-04
## 22310 HDI out VNM -0.0039676458 0.0028026661 1.164980e-03
## 2511 politicalstability out BLZ -0.0005327113 -0.0001304474 6.631587e-04
## 17711 politicalstability out SAU -0.0001422345 -0.0004458212 5.880557e-04
## 22311 politicalstability out VNM -0.0003940524 0.0001300490 2.640034e-04
## 2512 percent_foreign out BLZ -0.0129961436 0.0081264158 4.869728e-03
## 17712 percent_foreign out SAU -0.0038130784 0.0030500834 7.629950e-04
## 22312 percent_foreign out VNM -0.0107688054 0.0073338855 3.434920e-03
## 2513 percent_outside out BLZ 0.0186482168 -0.0121699888 -6.478228e-03
## 17713 percent_outside out SAU 0.0056955467 -0.0051847054 -5.108413e-04
## 22313 percent_outside out VNM 0.0144834736 -0.0100096678 -4.473806e-03res_impacts %>%
filter(model == "in",
iso3 %in% c("LBN", "CUB", "TCD"))## var model iso3 se1 se2 se3
## 47 hwf in CUB -0.0008519243 0.0010222358 -1.703115e-04
## 114 hwf in LBN -0.0006038206 0.0007015303 -9.770971e-05
## 200 hwf in TCD -0.0009255655 0.0010951575 -1.695920e-04
## 471 dry in CUB 0.0003217134 -0.0003825040 6.079058e-05
## 1141 dry in LBN 0.0002322123 -0.0002616565 2.944422e-05
## 2001 dry in TCD 0.0003497091 -0.0004097254 6.001628e-05
## 472 Lpop in CUB -0.0012732833 0.0024056921 -1.132409e-03
## 1142 Lpop in LBN 0.0002874049 0.0016205189 -1.907924e-03
## 2002 Lpop in TCD -0.0013343847 0.0026219778 -1.287593e-03
## 473 HDI in CUB 0.0037207092 -0.0043719031 6.511939e-04
## 1143 HDI in LBN 0.0027353794 -0.0029659142 2.305349e-04
## 2003 HDI in TCD 0.0040441014 -0.0046785193 6.344178e-04
## 474 politicalstability in CUB 0.0018190360 -0.0023355425 5.165065e-04
## 1144 politicalstability in LBN 0.0010732216 -0.0015982545 5.250329e-04
## 2004 politicalstability in TCD 0.0019584725 -0.0024933422 5.348697e-04
## 475 percent_foreign in CUB 0.0134033395 -0.0163236962 2.920357e-03
## 1145 percent_foreign in LBN 0.0088728257 -0.0108477550 1.974929e-03
## 2005 percent_foreign in TCD 0.0144502308 -0.0173739489 2.923718e-03
## 476 percent_outside in CUB -0.0171061830 0.0205420059 -3.435823e-03
## 1146 percent_outside in LBN -0.0125766247 0.0146503728 -2.073748e-03
## 2006 percent_outside in TCD -0.0187106727 0.0221576828 -3.447010e-03Compute the Moran plot on the residuals
We assess the presence of spatial autocorrelation in the residuals by computing the Moran plot, as illustrated in Figure D.12 of the article. This diagnostic helps evaluate whether the linear model adequately accounts for spatial dependence.
W_listw <- nb2listw(nb_out)
mat_out <- listw2mat(nb2listw(nb_out))
pal_HH <- RColorBrewer::brewer.pal(4, "Set2")
for(type in c("out", "in")) {
#pdf(paste0("figures/moran_residuals_", type, "flows.pdf"), width = 10, height = 7)
par(mfrow = c(2, 2), oma = c(0, 0, 0, 0), las = 1,
mgp = c(2.2, 1, 0))
for(r in c("ilr_1", "ilr_2")) {
# component 1
my_res <- residuals(my_lm[[type]][[r]])
#my_res <- temp_lm$out_ilr_1
# hwf + dry + Lpop + HDI + politicalstability + class_conflict +
#new_lm <- lm(out_ilr_1 ~ percent_foreign, data = temp_lm)
#my_res <- residuals(new_lm)
#lm.morantest(new_lm, nb2listw(nb_out))
my_res_lag <- lag.listw(nb2listw(nb_out), my_res)
HH <- my_res > 0 & my_res_lag > mean(my_res_lag)
HL <- my_res > 0 & my_res_lag <= mean(my_res_lag)
LH <- my_res <= 0 & my_res_lag > mean(my_res_lag)
LL <- my_res <= 0 & my_res_lag <= mean(my_res_lag)
col_vec <- rep(pal_HH[2], 230)
col_vec[HL] <- pal_HH[1]
col_vec[LH] <- pal_HH[3]
col_vec[LL] <- pal_HH[4]
# significant
local_ilr_1 <- localmoran(my_res, listw = W_listw)
col_vec <- ifelse(local_ilr_1[, 5] < 0.05, col_vec, scales::alpha(col_vec, 0.2))
cat("p.value: ", lm.morantest(my_lm[[type]][[r]], nb2listw(nb_out))$p.value, "\n")
par(mar = c(3.5, 4, 0.5, 0.5))
plot(my_res, my_res_lag,
xlab = paste0("Residuals from the ", type, "flows OLS model (", r, ")"),
ylab = "W x residuals",
col = col_vec, pch = 16)
abline(h = mean(my_res_lag, na.rm = T), lty = 2, col = "lightgrey")
abline(v = 0, lty = 2, col = "lightgrey")
abline(coefficients(lm(my_res_lag ~ my_res)), col = "red", lty = 1)
text(my_res[local_ilr_1[, 5] < 0.05], my_res_lag[local_ilr_1[, 5] < 0.05],
world_iso_3[local_ilr_1[, 5] < 0.05, ]$"VISUALIZATION_NAME", pos = 2,
cex = 0.6)
par(mar = c(0, 0, 0, 0))
plot(st_geometry(world_iso_3), col = col_vec, border = "lightgrey")
}
#dev.off()
}## p.value: 0.325363## p.value: 3.608411e-07
## p.value: 0.0188535## p.value: 0.3808837
Spatial CoDa Model
We now turn to the Spatial Autoregressive (SAR) model, adapted to compositional dependent variables (CoDa). This specification explicitly accounts for spatial autocorrelation in migration flows and provides a more reliable framework than the simple linear model.
# results modeling
res_sar_N <- vector("list", 2)
names(res_sar_N) <- c("out", "in")
# explanaotory variables
Xe <- temp_lm[, c(var_climate, var_controls)]
Xe$class_conflict_medium <- ifelse(Xe$class_conflict == "medium", 1, 0)
Xe$class_conflict_high <- ifelse(Xe$class_conflict == "high", 1, 0)
Xe$class_conflict_very_high <- ifelse(Xe$class_conflict == "very high", 1, 0)
Xe$class_conflict <- NULL
Xe <- data.frame(constant = 1, Xe)
nom_vars <- names(Xe)
Xe <- as(Xe, "matrix")
# dependent variables
Ye <- vector("list", 2)
names(Ye) <- c("out", "in")
for(type in c("out", "in")) {
Ye[[type]] <- as(ilr(temp_lm[,
paste0(type, "_", c("outwards", "returns", "transit"))],
V = ilr_base), "matrix")
res_sar_N[[type]] <- estimate_spatial_multi_gen_N(Y = Ye[[type]],
X = as(Xe, "matrix"),
W = OW,
method = "s2sls",
ind_RHO = matrix(c(T, T, T, T), 2, 2), compute_sd = T)
}We use the function sp_semi_elast_Ycompo() to compute the semi-elasticities, focusing on the subset of selected countries. The results correspond to Figures 4.8 and 4.9 of the thesis.
require(cartography)
ind_cnty <- list(
"out" = c(ind_viet, ind_saudi, ind_belize),
"in" = c(ind_chad, ind_cuba, ind_lebanon)
)
names_cnty <- list(
"out" = c("Viet Nam", "Saudi Arabia", "Belize"),
"in" = c("Chad", "Cuba", "Lebanon")
)
my_beta <-vector("list", 2)
my_RHO <-vector("list", 2)
names(my_beta) <- names(my_RHO) <- c("out", "in")
for(type in c("out", "in")) {
my_beta[[type]] <- res_sar_N[[type]]$est_beta
my_RHO[[type]] <- res_sar_N[[type]]$est_RHO
for(index_b in 2:3) {
elast_sp <- sp_semi_elast_Ycompo(Xe, my_beta[[type]],
index_b = index_b, method = "exact",
W = OW, RHO = my_RHO[[type]], V_y = ilr_base)
# We represent the semi-elasticities in the three countries chosen in the article
# fix the bounds of the colors
elas_k <- c(elast_sp[ind_cnty[[type]], , 1],
elast_sp[ind_cnty[[type]], , 2],
elast_sp[ind_cnty[[type]], , 3])
pos_class <- classInt::classIntervals(elas_k[elas_k>0], 6, "kmeans")
neg_class <- classInt::classIntervals(elas_k[elas_k<0], 6, "kmeans")
# we plot the se on the map
rwb <- colorRampPalette(colors = c("blue", "white", "red"))
#pdf(paste0("figures/elasticites_spatial_", type, "flows_" ,
# nom_vars[index_b], ".pdf"),
# width = 10, height = 7.5)
op <- par(mfrow = c(3, 3), oma = c(0, 1, 0, 0), mar = c(2, 2, 2, 2),
mgp = c(0, 0, 0))
for(l in 1:3) {
k <- ind_cnty[[type]][l]
world_iso_3$elas_1 <- elast_sp[k, , 1]
world_iso_3$elas_2 <- elast_sp[k, , 2]
world_iso_3$elas_3 <- elast_sp[k, , 3]
# we add one class for positive and negative
plot(st_geometry(st_buffer(world_iso_3[k, ], 1500000)), lty = 0,
ylab = paste0("Impacts due to a change \n of ",
nom_vars[index_b]," in ",
names_cnty[[type]][l]))
plot(sea, col = scales::alpha("lightblue", 0.1),
border = rgb(0.15, 0.15, 0.15), add = T)
choroLayer(world_iso_3, var = "elas_1",
breaks = c(neg_class$brks[1:6], pos_class$brks[2:7]),
col = rwb(11), legend.pos = "n", legend.values.rnd = 6,
legend.title.txt = "se (outward)", border = rgb(0.5, 0.5, 0.5),
lwd = 0.2, add = T)
if (l == 1) {
title("se Outward")
}
plot(st_geometry(st_buffer(world_iso_3[k, ], 1500000)), lty = 0)
plot(sea, col = scales::alpha("lightblue", 0.1),
border = rgb(0.15, 0.15, 0.15), add = T)
choroLayer(world_iso_3, var = "elas_2",
breaks = c(neg_class$brks[1:6], pos_class$brks[2:7]),
col = rwb(11), legend.pos = "n", legend.values.rnd = 6,
legend.title.txt = "se (returns)", border = rgb(0.5, 0.5, 0.5),
lwd = 0.2, add = T)
if (l == 1)
title("se Returns")
plot(st_geometry(st_buffer(world_iso_3[k, ], 1500000)), lty = 0)
plot(sea, col = scales::alpha("lightblue", 0.1),
border = rgb(0.15, 0.15, 0.15), add = T)
choroLayer(world_iso_3, var = "elas_3",
breaks = c(neg_class$brks[1:6], pos_class$brks[2:7]),
col = rwb(11), legend.pos = "right",
legend.values.rnd = 6, legend.title.txt = "semi-elasticities",
border = rgb(0.5, 0.5, 0.5), lwd = 0.2, legend.values.cex = 0.8,
add = T)
if (l == 1)
title("se Transit")
}
par(op)
#dev.off()
}
}## Direct impact:
## 0.002306955 -0.003754366 -0.003571608
## Indirect impact:
## -0.0001594533 0.0002624063 0.0002419488
## Total impact:
## 0.002147502 -0.00349196 -0.003329659
## Direct impact:
## -0.001240408 0.002056009 0.001855225
## Indirect impact:
## 7.639511e-05 -9.744572e-05 -0.0001653729
## Total impact:
## -0.001164013 0.001958564 0.001689852
## Direct impact:
## -0.001658768 0.002710886 -0.00186506
## Indirect impact:
## 0.0002477139 -0.000277907 -0.0001363335
## Total impact:
## -0.001411054 0.002432979 -0.002001393
## Direct impact:
## 0.0009010001 -0.001315185 0.0004995201
## Indirect impact:
## -1.750973e-05 7.942579e-05 -0.0001860443
## Total impact:
## 0.0008834903 -0.001235759 0.0003134758
Grouping the semi-elasticities
We use the function map_elas_ternary() to compute and
visualize the grouped semi-elasticities.
Figure 4.10 shows the direct (DE, red) and indirect (IE, blue)
semi-elasticities of migration compositions with respect to HWF. The top
panels display the results for outflows and the bottom panels for
inflows, with ternary diagrams on the left and barplots of aggregated
effects on the right.
coords_s <- st_coordinates(st_centroid(world_iso_3))
rownames(coords_s) <- world_iso_3$iso3#pdf("figures/arrow.pdf", width = 7, height = 6)
par(mfrow = c(2, 2), oma = c(0, 1, 0, 1), mar = c(0, 0.7, 0.1, 0.0))
for(type in c("out", "in")) {
for(index_b in c(2, 3)) {
map_elas_ternary(coords_s, add = F, row_to_plot = ind_cnty[[type]][1],
X = Xe, b_star = my_beta[[type]], index_b = index_b, W = OW,
type = 1, scale_ternary = TRUE,
RHO = my_RHO[[type]], delta = 50,
ternary_unit = T, range = 1,
cex.pts = 0.2, length = 0.05,
cex.legend = 0.6, lwd = 1.6,
legend.ternary = c("Outward", "Return", "Transit"),
names_ind = names_cnty[[type]][1], V_y = ilr_base)
map_elas_ternary(coords_s, add = T, row_to_plot = ind_cnty[[type]][2],
X = Xe, b_star = my_beta[[type]], index_b = index_b, W = OW,
type = 1, scale_ternary = FALSE,
RHO = my_RHO[[type]], delta = 50,
ternary_unit = T, range = 1,
cex.pts = 0.2, length = 0.05,
cex.legend = 0.6, lwd = 1.6,
legend.ternary = c("", "", ""),
names_ind = names_cnty[[type]][2], V_y = ilr_base)
map_elas_ternary(coords_s, add = T, row_to_plot = ind_cnty[[type]][3],
X = Xe, b_star = my_beta[[type]], index_b = index_b, W = OW,
type = 1, scale_ternary = FALSE,
RHO = my_RHO[[type]], delta = 50,
ternary_unit = T, range = 1,
cex.pts = 0.2, length = 0.05,
cex.legend = 0.6, lwd = 1.6,
legend.ternary = c("", "", ""),
names_ind = names_cnty[[type]][3], V_y = ilr_base)
if(type == "out" & index_b == 2) {
mtext("Outflows model", side = 2)
title("HWF ",
line = -1, cex.main = 1, font.main = 2)
}
if(type == "out" & index_b == 3) {
title("Dry ",
line = -1, cex.main = 1, font.main = 2)
}
if(type == "in" & index_b == 2)
mtext("Inflows model", side = 2)
}
}
legend("topright", legend = c("direct", "indirect"), pch = c(NA, NA),
lwd = 3, col = c("#E16A86", "#009ADE"),
cex = 0.9, lty = NA, box.lwd = -1) #normal legend. Do not plot line
par(font = 5) #change font to get arrows
legend("topright", legend = c(" ", " "),
pch = c(174, 174),
lwd = 3, col = c("#E16A86", "#009ADE"),
cex = 0.9, lty = NA, bty = "n", box.lwd = -1) 
#dev.off()We summarize the results into direct, indirect and total effects:
spatial_direct <- apply(elast_sp, 3, function(x) diag(x))
spatial_total <- apply(elast_sp, c(1, 3), sum)
spatial_indirect <- spatial_total - spatial_direct
colnames(spatial_direct) <- colnames(spatial_total) <-
colnames(spatial_indirect) <- c("Outward", "Return", "Transit")We now plot the individual barplot
source("codes/map_elas_bar.R")
#pdf("figures/bar_local_cantons.pdf",
# width = 16 * 0.5, height = 10 * 0.45)
par(oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0),
mfrow = c(2, 2))
for(type in c("out", "in")) {
for(index_b in c(2, 3)) {
elast_sp <- sp_semi_elast_Ycompo(Xe, my_beta[[type]],
index_b = index_b, method = "appro",
delta = 0.001,
W = OW, RHO = my_RHO[[type]], V_y = ilr_base)
spatial_direct <- apply(elast_sp, 3, function(x) diag(x))
spatial_total <- apply(elast_sp, c(1, 3), sum)
spatial_indirect <- spatial_total - spatial_direct
colnames(spatial_direct) <- colnames(spatial_total) <-
colnames(spatial_indirect) <- c("Outward", "Return", "Transit")
plot(sea, col = scales::alpha("lightblue", 0.2),
border = rgb(0.15, 0.15, 0.15))
if(type == "out" & index_b == 2)
title("HWF", line = -1.1)
if(type == "out" & index_b == 3)
title("Dry", line = -1.1)
plot(world_iso_3[ind_cnty[[type]], ], add = T,
border = "red", lwd = 1, lty = 2, col = "yellow")
map_elas_bar(spatial_direct, spatial_indirect, spatial_total,
coords_s, add = T,
index_to_plot = ind_cnty[[type]], plot_scale = F,
q4 = pal_col_ternary[c(1, 5, 9)],
width_bar = 1.3, max_bar = 5, cex.legend = 0.5, x.legend = "topleft",
label_s = F, round.scale = 0,
print_legend = ifelse(type == "out" & index_b == 2, T, F))
#text(coords_s[ind_cnty[[type]], 1] - 1500000, coords_s[ind_cnty[[type]], 2],
# world_iso_3[ind_cnty[[type]], ]$VISUALIZATION_NAME,
# pos = 3, cex = 0.6, col = rgb(0.1, 0.1, 0.1))
if(type == "out" & index_b == 2)
mtext(side = 2, text = "Outflows model", line = -1.1)
if(type == "in" & index_b == 2)
mtext(side = 2, text = "Inflows model", line = -1.1)
}
}## Direct impact:
## 0.002306956 -0.00375436 -0.003571603
## Indirect impact:
## -0.0001594533 0.0002624064 0.0002419488
## Total impact:
## 0.002147503 -0.003491954 -0.003329654## Direct impact:
## -0.001240408 0.002056011 0.001855226
## Indirect impact:
## 7.639511e-05 -9.744571e-05 -0.0001653729
## Total impact:
## -0.001164013 0.001958565 0.001689853## Direct impact:
## -0.001658768 0.002710888 -0.001865059
## Indirect impact:
## 0.0002477139 -0.000277907 -0.0001363335
## Total impact:
## -0.001411054 0.002432981 -0.002001393## Direct impact:
## 0.0009010002 -0.001315184 0.00049952
## Indirect impact:
## -1.750974e-05 7.942578e-05 -0.0001860443
## Total impact:
## 0.0008834905 -0.001235758 0.0003134757
#dev.off()index_b <- 2
my_delta <- 50
# merge spatial data
for(type in c("out", "in")) {
#pdf(paste0("figures/hwf_semi_", type , ".pdf"),
# width = 13 * 0.8, height = 4.35 * 0.8)
par(oma = c(0, 0, 0, 0), mfrow = c(1, 2))
# ternary
par(mar = c(0.2, 0, 0.2, 0))
map_elas_ternary(coords_s, add = F, row_to_plot = ind_cnty[[type]][1],
X = Xe, b_star = my_beta[[type]], index_b = index_b, W = OW,
type = 1, scale_ternary = TRUE,
RHO = my_RHO[[type]], delta = my_delta,
ternary_unit = T, range = 1,
cex.pts = 0.2, length = 0.05,
cex.legend = 0.8, lwd = 1.6,
legend.ternary = c("Outward", "Return", "Transit"),
names_ind = names_cnty[[type]][1], V_y = ilr_base)
map_elas_ternary(coords_s, add = T, row_to_plot = ind_cnty[[type]][2],
X = Xe, b_star = my_beta[[type]], index_b = index_b, W = OW,
type = 1, scale_ternary = FALSE,
RHO = my_RHO[[type]], delta = my_delta,
ternary_unit = T, range = 1,
cex.pts = 0.2, length = 0.05,
cex.legend = 0.6, lwd = 1.6,
legend.ternary = c("", "", ""),
names_ind = names_cnty[[type]][2], V_y = ilr_base)
map_elas_ternary(coords_s, add = T, row_to_plot = ind_cnty[[type]][3],
X = Xe, b_star = my_beta[[type]], index_b = index_b, W = OW,
type = 1, scale_ternary = FALSE,
RHO = my_RHO[[type]], delta = my_delta,
ternary_unit = T, range = 1,
cex.pts = 0.2, length = 0.05,
cex.legend = 0.6, lwd = 1.6,
legend.ternary = c("", "", ""),
names_ind = names_cnty[[type]][3], V_y = ilr_base)
if(type == "out")
mtext(side = 2, text = "Outflows model", line = -1.1)
if(type == "in")
mtext(side = 2, text = "Inflows model", line = -1.1)
# Map
#####################################.
par(mar = c(0, 0, 0, 0))
elast_sp <- sp_semi_elast_Ycompo(Xe, my_beta[[type]],
index_b = index_b, method = "appro",
delta = 0.001,
W = OW, RHO = my_RHO[[type]], V_y = ilr_base)
spatial_direct <- apply(elast_sp, 3, function(x) diag(x))
spatial_total <- apply(elast_sp, c(1, 3), sum)
spatial_indirect <- spatial_total - spatial_direct
colnames(spatial_direct) <- colnames(spatial_total) <-
colnames(spatial_indirect) <- c("Outward", "Return", "Transit")
plot(sea, col = scales::alpha("lightblue", 0.2),
border = rgb(0.15, 0.15, 0.15))
plot(world_iso_3[ind_cnty[[type]], ], add = T,
border = "red", lwd = 1, lty = 2, col = "yellow")
map_elas_bar(spatial_direct, spatial_indirect, spatial_total,
coords_s, add = T,
index_to_plot = ind_cnty[[type]], plot_scale = F,
q4 = pal_col_ternary[c(1, 5, 9)],
width_bar = 1.3, max_bar = 5, cex.legend = 0.5, x.legend = "topleft",
label_s = F, round.scale = 0,
print_legend = T)
#dev.off()
}## Direct impact:
## 0.002306956 -0.00375436 -0.003571603
## Indirect impact:
## -0.0001594533 0.0002624064 0.0002419488
## Total impact:
## 0.002147503 -0.003491954 -0.003329654
## Direct impact:
## -0.001658768 0.002710888 -0.001865059
## Indirect impact:
## 0.0002477139 -0.000277907 -0.0001363335
## Total impact:
## -0.001411054 0.002432981 -0.002001393